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I.  INTRODUCTION 


In  a  recent  report  we  discussed  calculations  performed  on  a  model 
for  simulating  shock  propagation  in  a  one-dimensional,  discrete  lattice. 

The  technique  employed  was  to  solve  numerically  the  classical  equations 
of  motion  of  the  atoms  in  the  chain  as  they  responded  to  the  shock  wave. 

From  the  solutions  of  the  equations,  one  can  then  infer  the  macroscopic 
properties  of  the  shock-compressed  lattice  by  taking  appropriate  aver- 

2 

ages.  This  technique^  is  referred  to  as  computer  molecular  dynamics 

anil  it  is  one  of  the  few  feasible  methods  for  studying  nonlinear,  non-  i 

equilibrium,  many-body  problems. 

In  Reference  1  we  emphasized  that  the  one-dimensional  calculations 
represented  only  an  initial  effort  and  that  we  intended  to  develop  a 

computer  program  for  solving  the  full  three-dimensional  problem.  The  , 

present  communication  is  the  first  in  a  series  of  two  reports  in  which 

we  will  discuss  the  three-dimensional  calculation.  In  the  current  work,  / 

we  discuss  the  model  very  briefly  and  confine  most  of  our  attention  to 

the  significant  results  of  the  calculation.  In  the  second  report'1, 
the  deta i 1 s  of  the  calculation,  including  a  discussion  of  the  computer 
program,  are  presented. 

The  need  for  a  microscopic,  discrete- lattice  treatment  of  shock 
propagation  and  its  relevance  to  Army-related  problems  have  been  thor¬ 
oughly  discussed  in  Reference  1  and  will  not  be  discussed  here.  Before 
proceeding,  however,  it  is  perhaps  worthwhile  to  indicate  briefly  some 
background  on  the  three-dimensional  problem  and  the  rationale  for  our 
specific  calculations.  Three-dimensional  CMI)  studies  of  shock  propa- 
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gat  ion  in  crystal  lattices  have  been  undertaken  by  Tsai  and  coworkers 


1.  .I.D.  Powell  and  .1.11.  Bat t eh,  "Shock  Propagation  in  the  One- 
Dimensional  Lattice",  BR1,  Report  No.  2009,  US  Army  Ballistic 
Research  Laboratories,  Aberdeen  Proving  Ground,  Ml),  1977.  See 
also  J.ll.  Batteh  and  .J.D.  Powell,  "Shock  Propagation  in  the  One- 
Dimensional  Lattice  at  Nonzero  Initial  Temperature",  .J.  Appl. 

Phys.  £9,  3933  (1978)  and  .J.ll.  Batteh  and  .J.D.  Powell,  in 
Soli  tons  in  Action,  edited  by  K.  Lonngren  and  A.  Scott  (Academic 
Press,  New  York,  1978,  p.  257. 

2.  B..J.  Alder  and  T.E.  Wainwright,  "Studies  in  Molecular  Dynamics. 

I.  (leneral  Method",  ,J.  Chem.  Phys.  3J_,  459  (  1959). 

3.  .J.D.  Powell  and  .J.ll.  Batteh,  "Shock  Propagation  in  the  Three- 
Dimensional  Lattice.  II.  Method  of  Calculation",  BRL  Report 
(to  be  published  concurrently  with  this  report). 

4.  R.A.  MacDonald  and  D.ll.  I'sai,  "Molecular  Dynamical 
of  Energy  Transport  in  Crystalline  Solids",  Phys. 

(1978). 

5.  D.ll.  Tsai  and  C.W.  Beckett,  "Shock  Wave  Propagation 
Lattices",  .J .  Ccophys.  Res.  7^,  2601  (1966). 

6.  D.ll.  Tsai  and  R.A.  MacDonald,  "Second  Sound  in  a  Solid  Under  Shock 
Compression",  .) .  Phys.  C  6,  1.171  (1973). 
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over  the  last  several  years.  An  important  and  controversial  conclusion 
of  their  work  has  been  that  the  shock  profile  is  not  steady  in  time  as 
is  generally  assumed  in  continuum  calculations.  Rather  than  observing 
a  steady  state,  these  workers  concluded  that  the  edge  of  the  thermally 
equilibrated  region  behind  the  shock  front  propagates  more  slowly  into 
the  lattice  than  does  the  front  itself.  Consequently,  the  transition 
region  between  the  front  and  the  region  of  equilibrium  grows  with  time. 
This  anomalous  behavior  was  explained  by  Tsai  as  being  a  natural 

7  8 

extension  of  the  low-temperature  phenomenon  of  second  sound  ’  .  For 
reasons  discussed  in  Reference  1,  however,  this  explanation  is  not 
entirely  satisfactory. 

In  addition  to  the  work  of  Tsai,  similar  calculations  have  been 

9-11 

performed  by  Paskin  and  coworkers  .  Although  both  investigators 
appear  to  obtain  similar  numerical  data,  substantial  disagreement  exists 
regarding  the  interpretation.  In  particular,  Paskin  has  disputed  Tsai's 
nonsteady  interpretation;  in  fact,  he  claims  essential  agreement  with 
the  results  of  continuum  models. 

12-19 

In  Reference  1  and  in  other  calculations  ,  both  analytic  and 
numerical,  shock  propagation  in  a  one-dimensional  lattice  has  been 


7.  J.C.  Ward  and  J.  Wilks,  "Second  Sound  and  the  Thermo-Mechanical 
Effect  at  Very  Low  Temperatures",  Phil.  Mag.  43,  48  (1952). 

8.  M.  Chester,  "Second  Sound  in  Solids",  Phys.  Rev.  13d,  2013  (1963). 

9.  A.  Paskin  and  G.J.  Dienes,  "Molecular  Dynamic  Simulations  of  Shock 
Waves  in  a  Three-Dimensional  Solid",  J.  Appl.  Phys.  4_3,  1605  (1972). 

10.  A.  Paskin  and  G.J.  Dienes,  "A  Model  for  Shock  Waves  in  Solids  and 
Evidence  for  a  Thermal  Catastrophe",  Solid  State  Comm.  Y7_,  19? 

(1975) . 

11.  A.  Paskin,  A.  Gohar,  and  G.J.  Dienes,  "Simulations  of  Shock  Waves 
in  Solids",  J.  Phys.  C  1£,  L563  (1977). 

12.  J.  Tasi,  "Perturbation  Solution  for  Growth  of  Nonlinear  Shock 
Waves  in  a  Lattice",  J.  Appl.  Phys.  4_3,  4016  (1972).  See  also 
Erratum  [J.  Appl.  Phys.  44,  1414  (1973)1- 

13.  J.  Tasi,  "Perturbation  Solution  for  Shock  Waves  in  a  Dissipative 
Lattice",  J.  Appl.  Phys.  44_,  2245  (1973). 

14.  J.  Tasi,  "Far-Field  Analysis  of  Nonlinear  Shock  Waves  in  a  Lattice", 
J.  Appl.  Phys.  4£,  4569  (1973). 

15.  G.F..  Duvall,  R.  Manvi,  and  S.C.  Lowell,  "Steady  Shock  Profile  in 
a  One-Dimensional  Lattice",  .J.  Appl.  Phys.  40,  3771  (1969). 

16.  R.  Manvi,  G.F..  Duvall,  and  S.C.  Lowell,  "Finite  Amplitude  Longi¬ 
tudinal  Waves  in  Lattices",  Int.  J.  Mech.  Sci.  U^,  1  (1969). 

17.  R.  Manvi  and  G.E.  Duvall,  "Shock  Waves  in  a  One-Dimensional,  Non- 
Dissipating  Lattice",  Brit.  J.  Appl.  Phys.  2,  1389  (1969). 

18.  .J.R.  Hardy  and  A.M.  Karo,  "Theoretical  Studies  of  Soliton-Like 
Behavior  of  Shocks  in  One-Dimensional  Systems",  LLL  Report  No. 
UCRL-79259,  Lawrence  Livermore  Laboratory,  Livermore,  CA,  1977. 

19.  B.L.  Holian  and  G.K.  Straub,  "Molecular  Dynamics  of  Shock  Waves 
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.will  res  soil .  1'hc  re  is  no  doubt  that  the  profile  is  nonsteady  in  one 

I  1  1 4 

dimension  and  some  investigators  ’  "  have  attributed  the  nonsteady 
behavior  to  the  propagation  behind  the  front  of  well-defined  pulses 
called  solitary  waves.  These  pulses  have  been  found  to  propagate  at 
speeds  which  vary  with  amplitude  and,  thus,  they  tend  to  spread  apart 
as  they  propagate.  The  spreading  effect,  then,  gives  rise  to  the 
nonsteady  behavior  observed  in  the  calculations.  The  stability  of  the 

pulses  was  also  studied*  and  they  were  found  to  be  stable  to  within 
the  numerical  accuracy  of  the  computer  data.  These  stable  pulses  are 

’() 

referred  to  as  soli  tons 

The  extent  to  which  solitary-wave  propagation  will  significantly 
affect  shock  propagation  (or  other  phenomena)  in  real  three-dimen¬ 
sional  crystals  is  not  currently  known.  It  is  conceivable,  for  instance, 
that  these  well-defined  pulses  will  not  remain  stable  when  they  en¬ 
counter  oscillations  transverse  to  their  direction  of  propagation.  In 
such  a  case  one  might  expect  that  their  influence  upon  the  shock  pro¬ 
file  would,  at  most,  be  a  transient  effect  in  a  three-dimensional 
model.  Because  of  the  importance  of  shock  propagation  to  such  Army- 
related  problems  as  detonation,  however,  it  is  of  obvious  importance 
to  settle  the  question  conclusively  as  well  as  to  check  the  general 
validity  of  the  continuum  approximations. 

The  purpose  of  the  calculations  carried  out  in  this  and  the  suc¬ 
ceeding  report  is,  therefore,  twofold:  first,  we  wish  to  determine  the 
extent  to  which  shock  propagation  in  a  three-dimensional  discrete 
latt  ice  is  at  variance  with  assumptions  made  in  continuum  models.  Of 
particular  interest  is  to  determine  whether  solitary  waves  propagate 
in  the  lattice  and,  if  so,  their  effects  upon  the  shock  profile. 

Second,  we  wish  to  develop  an  in-house  capability  for  performing  Oil) 
calculations  so  that  the  technique  can  be  applied  to  other  problems  of 
\rm\  interest. 

The  organization  of  this  first  report  is  as  follows.  In  Section  2 
we  describe  the  model  under  study,  write  down  the  equations  of  motion 
to  be  solved,  and  very  briefly  describe  the  method  of  solution.  Details 
of  the  calculation  arc  contained  in  Reference  3  which  is  devoted  to  that 
subject.  In  Section  3  the  numerical  results  of  the  calculation  are 
presented,  and  in  Section  I  the  shock  profiles  discussed,  finally, 
in  Section  5  we  discuss  the  results  in  greater  detail  and  compare  our 
calculations  with  those  of  other  investigators.  Some  indication  of 
possible  future  work  is  also  given. 

20.  N..I.  Zabusky  and  M.O.  kruskal,  "Interaction  of  Solitons  in  a 

Col  I i s ion  less  Plasma  and  the  Recurrence  of  Initial  States",  Pins. 

Rev.  Letters  15,  210  (IPOS). 


p 


2.  MODEL 


The  model  whose  properties  we  wish  to  study  consists  of  a  pure, 
monatomic,  FCC  lattice  which  is  made  as  long  as  necessary  in  the  z 
direction  and  which  is  periodic  in  the  x  and  y  directions.  The  peri¬ 
odicity  in  these  directions  is  characterized  by  the  integers  L  and  L  , 

x  y 

respectively.  Thus,  for  any  function  F  which  depends  upon  the  veloci¬ 
ties  and  displacements  of  the  atoms  in  the  lattice  we  have 


F(x+£LxaQ,  y+mlyxo,z)  =  F(x,y,z) 


(2.1) 


where  l  and  m  are  arbitrary  integers  and  a is  the  lattice  constant  or 
cube  edge  of  the  conventional  cell. 


The  atoms  will_be  assumed  to  interact  via  a  Morse-type  interatomic 
potential  and  therefore  the  Hamiltonian  of  the  lattice  can  be  written 

-R (A  |ft.  ..  -ft.  1-1)  "I  2 
j  -  o'  ljk  2,mn  1  i 

V  .  +  1/2  l  |e  -1  I  •  (2.2) 

i,j,k  J  i,j,k 

Si  ,m,n 


II  =  1/2  l 


[■ 


In  Eq.  (2.2)  all  quantities  have  been  made  nondimensional .  H  represents, 
of  course,  the  total  energy  in  the  lattice  and  has  been  normalized  by  D, 

the  dissociation  energy  of  a  single,  isolated  atom  pair;  V. is  the 

velocity  of  the  (ijk)  ‘  atom,  normalized  by  /n/m,  where  m  is  the  atomic 


mass 

V 


th 


;  R. is  the  position  vector  of  the  (ijk)  atom,  normalized  by 
ijk 

A  is  the  lattice  constant,  normalized  by  r  ,  the  separation  of  an 


isolated  atom  pair  at  minimum  potential;  and  R  is  a  dimensionless 
parameter  representing  the  degree  of  nonlinearity  in  the  Morse  potential. 
The  sum  over  (ijk)  runs  over  all  atoms  in  the  lattice  and  that  over 


th 


for  which 


(imn)  is  taken  over  all  atoms  in  the  vicinity  of  the  (ijk) 
the  potential  interaction  is  appreciable. 

The  equation  of  motion  satisfied  by  the  ( i j k ) atom  can  be  found 
in  a  straightforward  manner  from  Eq.  (2.2)  and  the  result  is 


ft.  ..  =  2  R  A  ft.  .. 
ijk  o  ijk 


(2.3) 


where  ft. is  the  nondimensional  force  (normalized  by  2  RD/r  )  exerted 

ljk  th  0 

on  the  (ijk)  atom  by  the  remaining  atoms  in  the  lattice  and  each  dot 

represents  differentiation  with  respect  to  the  dimensionless  time  t, 
i.e.,  the  real  time  normalized  by  /m/D  aQ. 
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Explicitly,  F. is  given  by 
1  j  R 


v  ■  ,  i 

£,m,n 


-2R(A  |ft.  ..40  1-1)  -R(A  j  R.  . ,  -R.  1-1) 

o  1  ilk  8-mn1  o'  ijk  limn  1 

e  J  -  e  J 


t. 

x  ijk  ftmn 


IV* 


r  2 . 4 ) 


X.mn 1 


In  all  our  calculations  we  will  be  concerned  with  solving  Eqs. 

(2.3)  numerically  to  determine  the  temporal  evolution  of  the  position 
and  velocity  of  each  atom  in  the  lattice  subject  to  some  specific  set 
of  initial  and  boundary  conditions.  The  initial  conditions  of  each 
atom  in  the  lattice  are  specified  randomly,  consistent  with  the  assump¬ 
tion  that  prior  to  compression  the  lattice  is  in  thermal  equilibrium  at 

t  b 

temperature  T.  Thus,  the  probability  that  the  a  Cartesian  component  of 


V. ..  lies  between  V  and  V+dV  is  given  by 
i  j  k  & 

function  0 


■  W 


1/2  -yV  /. 


d  V 


where  y  is  the  ratio  of  the  dissociation 
viz. , 


Y 


the  Maxwellian  distribution 


(2.5) 

energy  to  the  thermal  energy, 

(2.6) 


and  where  kD  is  Boltzmann's  constant.  The  boundary  condition  applied 

D 

to  the  model  is  such  that  the  first  plane  of  atoms  normal  to  the  z  axis, 
located  at  z=0,  is  driven  to  the  right  at  a  constant  nondimensional 
velocity  U  .  This  compression  produces  a  shock  wave  propagating  in  the 

z  direction  and,  from  the  solution  of  Eqs.  (2.3),  we  infer  the  nature 
of  the  response  of  the  lattice  to  the  shock  wave. 

3.  NUMERICAL  RESULTS 


In  all  our  calculations  we  have  assumed,  unless  otherwise  noted, 
that  the  anharmonicity  factor  R  in  Eq.  (2.2)  was  given  by  6.29,  a  value 
which  leads  to  a  reasonable  representation  for  a  lattice  of  nickel 
2 1 

atoms  .  Furthermore,  we  have  assumed  that  only  atoms  which  were 

TT.  F.  Milstein,  "Applicability  of  Exponentially  Attractive  and 

Repulsive  Interatomic  Potential  Functions  in  the  Description  of 
Cubic  Crystals",  J.  Appl.  Phys.  44,  3825  (1973). 
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separated  by  a  distance  of  unity  or  less  (real  distance  normalized  by 
a^)  exerted  an  appreciable  force  on  one  another.  This  ass'imption  is 

equivalent  to  assuming  that,  in  the  equilibrium  lattice  at  zero  tem¬ 
perature,  only  an  atom's  first-  and  second-nearest  neighbors  contribute 
significantly  to  its  potential  interaction.  The  assumption  was  found  to 
be  reasonable  for  the  currently  used  value  of  R.  The  lattice  constant 
was  then  calculated  (see  Reference  3)  and  found  to  be  1.4034.  The 
parameter  y  was  chosen  to  be  16  which  corresponds  to  a  lattice  at  roughly 
150°K.  The  small  value  for  the  initial  temperature  was  chosen  so  as  to 
minimize  the  possibility  of  the  lattice  melting  upon  compression.  Some 
calculations  have  also  been  done  with  y=9,  however,  which  corresponds  to 
an  initial  temperature  near  room  temperature,  and  no  qualitative  differ¬ 
ences  with  results  for  the  lower-temperature  calculations  were  observed. 

To  assist  in  interpreting  the  results  reported  in  the  following 
sections,  we  have  provided  in  Table  I  a  list  of  symbols  with  their 
corresponding  normalizing  factors.  The  specific  numerical  values  for 
the  parameters  used  in  our  simulation  of  a  nickel  lattice  are  °lso 
included. 

TABLE  I.  Symbols  and  Normalizing  factors. 


Quantity 

Symbo 1 

Normalizing  Factor 
or 

Numerical  Value 

Velocity  of  particle  (i,j,k) 

V 

/n/m 

Nondimensional  lattice  constant 

A 

r 

0 

o 

Position  of  particle  (i,j,k) 

R.  .. 
ljk 

a 

0 

Potential  energy  of  particle  (i,j,k) 

4.  .. 

ij  k 

D 

liens  ity 

P 

po 

Force  exerted  on  particle  (i,j,k) 

F>ijk 

2RD/r 

o 

Time 

T 

/m/li  a 

0 

*+ 

,  3 

Stress  tensor 

0 

4D/«o 

Compression  velocity 

u 

p 

/f)/m 

Temperature 

6 

DA 

°  -20 

3.5  X  10  j 

Dissociation  energy 

D 

Mass  of  particle 

m 

9.7  X  10"26kg 

Equilibrium  spacing  for  isolated 
pa i r  of  atoms 

r 

2.53  X  l(fll)m 

Lattice  constant 

a 

3.52  X  10~10m 

Initial  density  ahead  of  shock 

o 

po 

n  o  Y  ”3 

9.2  X  10  m 

Nonlinearity  parameter 

R 

6.  29 

3.1  The  Initially  Quiescent  Lattice  and  the  Existence  of 

Solitary  Waves. 

The  simplest  instructive  calculation  that  can  be  carried  out  is 
one  for  which  the  initial  temperature  of  the  lattice  is  zero  (y-*00) . 

In  that  case  each  atom  in  the  lattice  is  initially  at  rest  in  its 
equilibrium  position  and  the  shock  wave  produces  motion  only  along  the 
z  direction.  The  effect  of  the  shock  wave  can  be  determined  most  easily 
by  plotting  the  velocity-time  trajectories  of  various  planes  of  atoms, 
normal  to  the  z  axis,  as  they  are  encountered  by  the  shock. 

A  set  of  such  trajectories  is  shown  in  Figure  1  for  the  case  U  =2.0. 

Each  plot  is  begun  at  time  t  when  the  shock  wave  first  excites  the 

plane  in  question  and  the  notation  V  denotes  the  (common]  velocity  of 

atoms  in  the  a  plane.  As  can  be  seen  from  the  figure  the  shock  causes 
the  second  plane  of  atoms  to  move  initially  along  an  oscillatory-wave 
profile  similar  to  that  characteristic  of  atoms  in  a  one-dimensional, 

harmonic  chain^.  As  the  wave  propagates  farther,  however,  and  the  atoms 
become  farther  displaced  from  their  equilibrium  positions,  nonlinear 
effects  become  increasingly  important.  The  nonlinearity  of  the  lattice 
tends  to  steepen  the  wave  profiles  until  such  time  as  the  nonlinear 
effect  is  exactly  balanced  by  the  dispersive  effect  of  the  lattice  which 
tends  to  broaden  the  profile.  At  this  time  rather  well-defined  pulses 
known  as  solitary  waves  propagate  in  the  vicinity  of  the  shock  front. 

The  evolution  of  the  solitary  waves  can  be  seen  in  the  velocity¬ 
time  trajectories  of  the  20^'  and  40t*1  planes  in  the  figure.  At  the 
ti  h 

20  plane,  the  amplitudes  of  the  pulses  are  seen  to  decrease  as  one 
moves  away  from  the  shock  front.  Since  the  pulse  propagation  velocity 
decreases  with  decreasing  amplitude,  a  spreading  effect  ensues.  Conse¬ 
quently,  at  the  40t*1  plane  the  pulses  are  seen  to  be  slightly  more 

distantly  spaced  than  when  the  shock  front  was  at  the  20  plane.  As 
the  pulses  near  the  front  propagate  still  farther  into  the  lattice,  it 
is  found  that  their  shapes  remain  essentially  constant  and  they  are 
then,  by  definition,  solitary  waves. 

In  Reference  1  we  have  shown  that  for  a  one-dimensional,  Morse- 
potential  lattice,  solitary  waves  propagate  near  the  shock  front  not 
only  for  the  initially  quiescent  lattice,  but  also  for  the  initially 
thermalized  lattice.  In  the  latter  case,  the  solitary  waves  have 
varying  amplitudes  near  the  shock  front  and  the  spreading  effect  alluded 
to  earlier  prevents  the  shock  profile  from  approaching  a  steady  state. 
Furthermore,  we  have  investigated  the  stability  of  the  solitary  waves 
in  one  dimension  and  found  that,  to  within  our  numerical  accuracy, 
they  are  stable  to  mutual  collisions  and  to  interaction  with  the  thermal 
background  of  the  lattice.  One  can  expect,  therefore,  that  the  stability 


Figure  1.  Evolution  of  solitary  waves  at  the  shock  front  in  the 
initially  quiescent  lattice. 
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of  the  solitary  waves  will  prevent  the  rapid  establishment  of  thermal 
equilibrium  behind  the  shock  front. 

It  remains  to  determine  whether  similar  effects  occur  when  the 
three-dimensional,  initially  thcrinalized  lattice  is  subjected  to  shock 
compression.  Preparatory  to  this  investigation,  we  have  undertaken  a 

22 

rather  thorough  study'"'  of  the  properties  of  solitary  waves  in  the  three- 
dimensional  lattice.  The  principal  conclusions  of  that  study,  which 
will  be  useful  in  interpreting  the  current  work,  are  as  follows. 

1.  Longitudinal  solitary  waves  in  a  three-dimensional 
lattice  are  not  so  stable  to  mutual  collisions  as  those  in  a  one¬ 
dimensional  chain,  and  the  degree  of  stability  decreases  with  increasing 
amplitude  of  the  colliding  pulses.  Despite  the  lack  of  complete  sta¬ 
bility,  however,  the  pulses  do  maintain  their  integrity  fairly  well. 

2.  The  solitary  waves  are  extremely  stable  to  small, 
longitudinal,  planar  oscillations. 

3.  Longitudinal  solitary  waves  are  not  stable  to  small, 
transverse,  planar  oscillations  whenever  the  amplitude  of  the  initial 
solitary  wave  is  sufficiently  large.  Even  in  this  case,  however,  most 
of  the  energy  associated  with  the  initial  wave  remains  localized  in  the 
form  of  coupled  longitudinal  and  transverse  solitary  waves  which  propa¬ 
gate  at  the  same  speed. 

4.  Solitary  waves  appear  to  be  more  stable  to  random  thermal 
oscillations  than  to  coherent  planar  oscillations. 

3.2  Shock  Propagation  in  the  Initially  Thermalized  Lattice. 

In  an  effort  to  determine  how  the  shock  profile  is  affected  by 
thermal  background  in  three  dimensions,  we  have  subjected  an  initially 
thermalized  lattice  to  shock  compression  in  the  same  manner  as  before. 

The  compression  velocity  was  again  given  by  U^=2.0  and  the  parameter  y 

was  chosen  to  be  lb.  The  periodicity  of  the  lattice  was  given  by  L  =L  =4 

which  corresponds  to  32  (unique)  atoms  in  each  plane  normal  to  the  z 
axis. 

Again,  it  is  instructive  to  plot  velocity-time  trajectories  of 
individual  planes,  where  now  we  define  the  velocity  of  a  plane  to  be 
the  average  of  the  velocities  of  the  atoms  contained  within  the  plane. 

figure  2  shows  such  a  trajectory  for  the  2()f^  plane  subsequent  to  its 
being  excited  by  the  shock.  From  the  plot  of  the  z  component,  t 

22.  .J  .11.  Batteh  and  .1.1).  Powell,  "Solitary-Wave  Propagation  in  the 
Three-Dimensional  Lattice",  Phys.  Rev.  B  (in  press). 
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we  can  see  that  solitary  waves  are  apparently  forming  in  the  vicinity 
of  the  shock  front  just  as  in  the  initially  quiescent  lattice.  Comparing 
this  figure  with  figure  1,  however,  we  observe  that  the  amplitudes  of 
the  pulses  arc  slightly  less  than  in  the  quiescent  lattice  and  the  pro¬ 
files  are  less  well  defined.  The  difference  arises,  no  doubt,  because 
in  the  thermal ized  lattice  the  solitary  waves  continuously  interact 
with  the  thermal  background  of  the  lattice.  It  can  also  be  observed 
that  the  shock  front  has  given  rise  to  some  fairly  large-amplitude 
planar  oscillations  in  the  transverse  directions. 

By  the  time  the  shock  front  has  reached  the  80^'  plane  in  the 
lattice,  the  trajectory  of  which  is  shown  in  Figure  3,  the  profiles  of 
the  pulses  near  the  front  have  become  significantly  more  well  defined. 

The  amplitude  of  the  leading  pulse,  however,  has  decreased  by  about 
25%.  It  is  conceivable  that  part  of  the  energy  associated  with  the 
original  longitudinal  solitary  wave  became  coupled  into  the  thermal 
background  as  this  pulse  propagated  farther  into  the  lattice.  Apparently, 
however,  a  substantial  part  of  the  energy  has  been  transferred  to  a 
rather  well-defined  planar  oscillation  in  the  x  direction  as  can  be 
seen  in  the  plot  of  (VX())  .  This  transverse  pulse  can  also  be  seen  in 

x 

the  plot  of  ,  but  its  amplitude  is  considerably  less  than  when 

u  x  t, 

the  front  is  near  the  8fT  plane.  Since  the  pulse  which  propagates  in 
the  x  direction  seems  to  be  exactly  in  phase  with  the  leading  longitu¬ 
dinal  solitary  wave,  and  travelling  at  the  same  propagation  velocity, 
we  are  led  to  conclude  that  the  two  pulses  are  precisely  the  coupled 

22 

solitary  waves  which  we  have  observed  in  our  previous  work-  .  These 
coupled  waves  were  found  to  arise  whenever  longitudinal  solitary  waves 
of  sufficiently  high  amplitude  encountered  small,  transverse  planar 
oscillations.  A  similar  pulse  can  be  seen  beneath  the  second  longitu¬ 
dinal  solitary  wave  in  Figure  3. 

It  appears  surprising  that  the  coupled  waves  appear  only  in  the 
x  and  z  directions  with  little  transfer  of  energy  into  the  y  direction. 
This  unusual  condition  apparently  arose  because  of  the  precise  micro¬ 
scopic  state  of  the  lattice  when  the  excitation  occurred.  At  other 
times  and  in  other  calculations  we  have  observed  significant  oscillations 
in  all  three  directions. 


We  have  allowed  the  shock  front  to  propagate  approximately  100 
planes  into  the  lattice  discussed  above.  Because  of  the  immense  amount 
of  computer  time  needed  to  solve  all  the  equations  of  motion  for  such 
a  large  lattice,  however,  we  have  also  repeated  the  calculation  using 


a  lattice  for  which  I.  =1.  =2. 

x  y 


This  lattice  contains  8  atoms  per  plane 


and,  because  of  the  smaller  cross  section,  we  were  able  to  allow  the 


shock  to  propagate  considerably  farther.  In  fact,  in  our  longest 
calculations,  the  shock  front  has  propagated  more  than  300  planes  into 


the  lattice. 


In  Figure  I  we  have  shown  the  velocity-time  trajectory  of  the  240^ 
plane  of  atoms  for  the  shock  calculation  in  the  smaller  lattice.  The 

results  are  somewhat  similar  to  those  observed  at  the  plane  in  the 

previous  calculation  except  that  the  longitudinal  wave  profiles  in  the 
vicinity  of  the  front  are  still  more  clearly  defined  and  some  of  the 
energy  has  been  transferred  to  the  y  direction.  In  fact,  the  leading 
longitudinal  solitary  wave  appears  to  be  accompanied  by  transverse 
solitary  waves  in  both  the  x  and  y  directions,  but  the  shapes  of  the 
profiles  are  rather  poorly  defined  because  of  interaction  with  the 
thermal  background.  In  order  to  verify  that  the  pulses  we  observe  are 
indeed  solitary  waves,  it  is  desirable  to  separate  them  from  the  thermal 
background.  The  separation  can  be  accomplished  most  readily  by  a  tech- 
nique  suggested  in  Reference  1.  It  was  observed  there  that  the  high- 
amplitude  solitary  waves  at  the  head  of  the  shock  front  would  propagate 
at  a  faster  velocity  than  the  rate  at  which  thermal  energy  would  diffuse 
into  the  lattice.  Therefore,  in  order  to  separate  the  two  effects,  we 
instantaneously  stopped  the  compression  when  the  shock  front  had  reached 

the  243r<^  plane  in  the  lattice.  The  atoms  beyond  the  243rt^  plane,  those 
in  the  uncompressed  lattice,  were  then  placed  at  rest  in  their  equili¬ 
brium  positions  so  that  the  hot,  shock-compressed  lattice  was  directly 
next  to  a  cold  lattice.  One  should  then  be  able  to  observe  the  well- 
defined,  high-amplitude  solitary  waves  propagating  into  the  cold  lattice, 
leaving  the  thermal  background  behind. 

The  results  of  the  calculation  are  shown  in  Figure  5.  We  have 

f"  h 

plotted  the  velocity-time  trajectory  of  the  260  plane  in  the  lattice 

after  the  compression  was  stopped  at  the  243r  plane.  Three  extremely 
well-defined  pulses  can  be  seen  to  propagate  exactly  in  phase  in  the 
three  Cartesian  directions. 

We  have  now  observed  solitary-wave  propagation  in  the  vicinity  of 
the  shock  front  in  a  number  of  cases.  To  summarize,  our  understanding 
of  precisely  how  this  phenomenon  occurs  is  as  follows.  Subjecting  the 
lattice  to  shock  compression  produces  a  spectrum  of  longitudinal  soli¬ 
tary  waves  which  propagate  into  the  lattice  and  interact  with  the  thermal 
background.  When  these  solitary  waves  are  subjected  to  transverse 
planar  oscillations  (however  small)  they  accentuate  those  oscillations 
and,  for  sufficiently  strong  longitudinal  waves,  transverse  solitary 
waves  are  generated.  The  transverse  waves  are  coupled  to  the  original 
wave  and  propagate  in  phase  with  it.  The  solitary  waves,  however,  are 
not  completely  stable  to  the  thermal  oscillations  in  the  lattice  and 
one  can  expect  that  they  will  eventually  decay  as  they  propagate  into 
the  lattice.  The  energy  which  evolves  as  a  result  of  the  decay  then 
increases  the  random  thermal  energy  (temperature)  behind  the  shock 
front.  In  the  following  section  we  will  discuss  the  effects  of  solitary 
waves  upon  the  shock  profiles. 
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Figure  4.  Velocity-time  tra 


4.  SHOCK  PROF  I  l.liS  AND  THERMODYNAMI C  QUANTITIES 


In  this  section  we  wish  to  examine  the  shock  profiles  predicted 
from  our  numerical  data.  We  will  be  particularly  interested  in  deter¬ 
mining  how  the  propagation  of  solitary  waves  in  the  lattice  affects 
these  profiles.  In  calculating  all  profiles,  we  have  averaged  the 
thermodynamic  variable  in  question  over  regions  of  the  lattice  which 
were  microscopically  large  but  macroscopically  small.  In  most  calcu¬ 
lations  the  regions  were  chosen  to  be  16  planes  in  length  (along  the 

z  axis)  and  contained,  therefore,  128  atoms  for  the  case  L  =L  =2. 

x  y 

4.1  Density,  Velocity,  and  Stress. 


The  simplest  thermodynamic  variables  to  calculate  and  interpret 
are  the  density,  average  velocity,  and  stress.  Their  calculation 
has  been  discussed  in  some  detail  in  Reference  3  and  here  we  merely 
quote  the  appropriate  formulas: 


(4.1) 


<  V  >  =  —  Tv... 

nR  R  ljk 


(4.2) 
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i  ,  m ,  n 
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i j  k  i.mn 
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(4.3) 


In  Eq.  (4.1)  p  represents  the  density,  AZ^.  the  (average)  final  separation 

along  the  z  axis  of  the  planes  which  bound  the  region  under  considera¬ 
tion,  and  AZ.  the  initial  separation.  The  density  has  been  normalized 

by  the  initial  density.  In  Eq.  (4.2),  <  V  >  is  the  average  velocity 
and  the  notation  }'  denotes  that  sum  was  only  over  atoms  in  region  R; 

n  is  the  number  of  atoms  in  that  region.  Finally  the  two  terms  on  the 

right-hand  side  of  liq .  (4.3)  represent  the  kinetic  and  potential  contri- 

4-¥ 

but  ions  to  the  pressure  or  stress  tensor  o  .respect i vely .  The  quantity 

F...  denotes  the  force  exerted  on  atom  (i,j,k)  by  atom  (fc,m,n).  The 

i  j  k  ,  Jo  inn  ^ 

stress  has  been  normalized  by  the  factor  4D /a'Q. 

In  Figure  6  we  have  plotted  the  density,  the  z  component  of  the 
average  velocity,  and  the  zz  element  of  the  stress  tensor  as  a  function 


of  position  along  the  z  axis.  The  values  were  obtained  for  the  case 

U  =2.0  for  L  =1.  =2  and  when  the  shock  was  at  approximately  the  32()t*1 
p  x  y 

plane.  Here  x  refers  to  the  time  at  which  the  shock  reached  that 
1  o 


plane.  In  order  to  plot  the  stress  on  the  same  vertical  axis  as  the 

other  two  variables  we  have  actually  plotted  a  / 1 0 . 

z  z 


As  can  be  seen  from  the  figure,  all  quantities  are  rather  uniform 
except  near  the  shock  front  where  a  rapid  variation  occurs  within  a 
distance  of  a  few  planes.  Thus,  the  profiles  appear  rather  similar  to 
those  expected  from  a  continuum-mechanical  analysis.  The  velocity, 
of  course,  approaches  the  value  2.0,  the  driving  velocity,  behind  the 
front.  The  density,  which  is  unity  ahead  of  the  front,  approaches  the 
value  1.2  behind  the  front,  corresponding  to  a  compression  of  20%.  The 
zz  element  of  the  stress  tensor  increases  from  a  value  of  about  0.5 
ahead  of  the  front  to  a  value  of  about  25  behind  the  front.  The  xx  and 
yy  elements,  not  shown  in  the  figure,  are  both  equal  to  about  20  behind 
the  front. 


We  have  drawn  a  number  of  other  profiles,  corresponding  to  those 
above,  but  when  the  shock  front  was  at  some  other  location  in  the 
lattice.  The  profiles  have  essentially  the  same  shapes  and  therefore 
appear  to  be  steady  in  time.  However,  all  three  thermodynamic  quanti¬ 
ties  are  fairly  insensitive  to  small  changes  in  the  state  of  the  lattice 
and  are  not  a  good  indication  of  whether  steady  state  exists. 

4.2  Temperature  Profiles  and  the  Question  of  Steady  State. 

The  thermal -energy  density,  defined  by 


(4.4) 


is  more  sensitive  to  small  changes  in  the  state  of  the  lattice  and  is, 
therefore,  more  appropriate  for  examining  the  existence  of  steady  state 
Physically,  G  represents  'hat  portion  of  the  kinetic  energy  in  region  R 
which  is  not  associated  with  translation  of  the  region  as  a  whole.  For 
the  special  case  in  which  R  is  in  local  thermodynamic  equilibrium,  it 
corresponds  to  the  temperature  (normalized  by  the  factor  D/k^)  in  the 

region.  As  a  matter  of  convenience,  however,  we  wi 1 1  still  refer  to  9 
as  calculated  by  Eq.  (4.4)  as  the  temperature  even  though  the  region 
for  which  it  is  calculated  may  not  be  in  equilibrium. 

Figure  7  shows  three  temperature  profiles  corresponding  to  times 
when  the  shock  front  was  at  planes  258,  291  and  3 23.  As  can  be  seen 
from  the  profiles,  the  temperature  near  the  front  substantially  over¬ 
shoots  its  final  constant  value.  The  overshoot  is  followed  by  a  region 
in  which  thermal  relaxation  occurs  and,  finally,  a  region  of  nearly 
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constant  temperature  (roughly,  0.17  or  three  times  ambient).  To 
understand  this  behavior  it  is  convenient  to  rewrite  Hq.  (4.4)  as 
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where  V 


P 

ijk 


denotes 


the  average  velocity  of  the  plane, 


normal  to  the  z 


axis,  which  contains  the  atom  (i,j,k).  According  to  Hq.  (4.5),  there 
are  two  contributions  to  the  temperature  in  the  region.  The  first 
contribution  arises  from  the  fluctuations  in  the  velocities  of  the  atoms 
about  the  mean  velocity  of  the  planes  in  which  they  are  located.  The 
second  term  arises  from  the  fluctuations  in  the  velocities  of  the  planes 
themselves  about  the  mean  velocity  in  the  region  as  a  whole.  As  we 
noted  in  Section  .3,  the  solitary  waves  propagating  just  behind  the 
shock  front  induce  large-amplitude  planar  oscillations  and  it  is  precisely 
the  contribution  of  these  oscillations  to  the  second  term  in  hq.  (4.5) 
which  leads  to  the  overshoot  in  the  temperature.  As  one  moves  away 
from  the  front  towards  the  driven  end  of  the  lattice,  the  planar  oscil¬ 
lations  decrease  in  amplitude  and,  consequently,  the  temperature  de¬ 
creases.  Finally,  far  behind  the  shock  front  the  planar  oscillations 
have  decayed,  each  plane  is  moving  essentially  at  the  compression 
velocity,  and  the  temperature  in  that  region  of  the  lattice  is  due 
primarily  to  the  random  motion  of  the  atoms  within  the  planes. 


It  is  useful  to  compare  these  temperature  profiles  to  those  ob¬ 
tained  for  an  initially  quiescent  lattice  also  compressed  at  the  rate 
11  =2.0.  We  have  plotted  two  profiles  in  Figure  8,  corresponding  to  times 

when  the  shock  front  was  at  planes  162  and  226.  The  profiles  in  Figure 
8  are  qualitatively  similar  to  the  temperature  profiles  in  the  transition 
region  in  Figure  7.  This  results  from  the  fact  that,  in  both  cases,  the 
temperature  is  dominated  by  the  contribution  arising  from  the  planar 
oscillations.  It  is  also  evident  from  Figure  8  that  the  shock  profile 
in  the  initially  quiescent  lattice  is  not  steady  in  time.  In  fact,  the 
transition  region  has  increased  by  about  60  planes  in  the  time  required 

for  the  shock  to  propagate  from  the  162nc*  plane  to  the  226ttl  plane. 

As  we  have  suggested  before,  the  growing  transition  region  results  from 
the  spreading  of  pulses  of  various  amplitude,  as  they  evolve  into 
solitary  waves. 


It  is  not  so  easy  to  ascertain  from  the  profiles  in  Figure  7  whether 
the  transition  region  is  growing  or  not.  Unfortunately,  thermal  back¬ 
ground  sufficiently  obscures  the  part  of  the  profile  resulting  from  the 
solitary  waves  that  a  spreading  effect,  if  it  exists,  cannot  be  unequi¬ 
vocally  observed.  We  have  plotted  and  examined  in  detail  many  profiles 
simila.  to  those  in  Figure  7  and  were  unable  to  come  to  any  conclusion. 


shoulu  emphasize  that  absence  of  steady  state  in  one  dimension 
and  in  the  three-dimensional  quiescent  lattice  does  not  necessarily 
imply  that  the  profile  will  be  unsteady  in  the  three-dimensional,  ther¬ 
mal  ized  lattice.  The  principal  difference  is  that,  in  the  first  two 
cases,  there  is  no  appreciable  decay  of  the  solitary  waves  as  they 
propagate  into  the  lattice  and  the  spreading  effect  must  persist.  In 
the  last  case,  however,  the  solitary  waves  at  the  front  are  not  completely 
stable.  Therefore,  as  they  give  up  their  energy  to  the  thermal  back¬ 
ground  and  their  propagation  velocity  decreases,  they  may  be  overtaken 
by  higher-amplitude  solitary  waves  which  were  initially  behind  them.  It 
is  then  conceivable,  if  not  likely,  that  such  a  profile  would  appear 
steady  in  a  global,  macroscopic  sense. 

As  we  have  pointed  out  in  Section  1,  Tsai  has  reported  that  the 
temperature  profiles  are  not  steady  in  time.  Our  results  are  not  in 
direct  conflict  with  this  allegation  but  we  should  point  out  that  the 
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ambient  temperature  employed  by  Tsai  was  extremely  low  (~38°K).  One 
would  clearly  expect  that,  for  such  a  low  thermal  background,  the  soli¬ 
tary  waves  would  decay  exceedingly  slowly.  Consequently,  the  spreading 
effect  and  thus  the  nonsteady  behavior,  would  be  expected  to  persist 
for  very  long  times  indeed,  perhaps  longer  than  the  duration  of  Tsai's 
ca lculat ions . 


In  addition  to  plotting  profiles  for  the  total  temperature,  we  have 
also  plotted  components  of  the  temperature  associated  with  each  Cartesian 
direction.  The  components  are  given  by  a  relation  analogous  to  Eq.  (4.4), 
namely 
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where  3  here  denotes  the  Cartesian  component  under  consideration.  The 
components  are  shown  in  figure  9  for  the  case  corresponding  to  the  last 
of  the  graphs  in  Figure  7.  It  is  interesting  to  note  that  overshoots 
occur  not  only  in  the  z  direction  but  also  in  the  x  and  y  directions. 
Temperature  overshoots  in  directions  normal  to  the  direction  of  shock 
propagation,  which  are  somewhat  unexpected,  result  from  the  transverse 
planar  oscillations  associated  with  the  coupled  solitary  waves  discussed 
in  Section  3.  It  is  interesting  that  we  did  not  observe  overshoots  in 
the  transverse  directions  for  more  weakly  shocked  lattices,  say,  for 
Up=1.0.  We  believe  that  their  absence  results  from  the  fact  that  the 

accompanying  longitudinal  solitary  waves  were  not  sufficiently  strong 

to  generate  the  transverse  waves.  In  fact,  we  have  demonstrated  in 
22 

previous  work  “  that  there  exists  a  threshold  amplitude  for  the  longi¬ 
tudinal  solitary  wave  below  which  coupled  solitary  waves  do  not  exist. 


4.3  Investigation  of  Equilibrium. 


In  the  preceding  section  we  have  used  the  term  "temperature"  rather 
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loosely  to  denote  the  nondirected  kinetic  energy  in  a  region  and  now  we 
wish  to  examine  the  extent  to  which  the  region  behind  the  shock  is  act¬ 
ually  in  thermal  equilibrium.  For  a  region  to  be  in  local  thermodynamic 
equilibrium,  it  is  necessary  that  (1)  the  thermal  energy  in  the  region 
be  equally  distributed  among  the  three  temperature  components,  (2)  the 
value  of  the  thermal  energy  differ  only  slightly  from  its  value  in 
adjacent  regions,  and  (3)  the  velocity  distribution  function  in  the 
region  be  Maxwellian.  Figure  9  clearly  demonstrates  that  the  first  two 
conditions  are  not  satisfied  just  behind  the  shock  front.  The  absence 
of  thermodynamic  equilibrium  in  this  region  is  not  unexpected  and  is  a 
direct  consequence  of  the  planar  oscillations  associated  with  the  propa¬ 
gation  of  the  solitary  waves.  On  the  other  hand,  far  behind  the  front 
where  the  planar  oscillations  have  decayed,  the  temperature  is  rela¬ 
tively  constant  and  the  thermal  energy  is  equally  partitioned,  suggesting 
that  there  the  lattice  may  be  in  equilibrium. 

In  an  effort  to  verify  these  observations,  we  have  plotted  a  number 
of  velocity  distribution  functions  for  various  regions  of  the  lattice 
at  several  times  during  the  shock-wave  calculations.  Typical  of  the 
distributions  obtained  are  those  shown  in  Figure  10,  corresponding  to 
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time  i  =12.80  when  the  shock  was  at  the  323  plane.  Plotted  on  the 
o  1 

vertical  axis  is  the  fractional  number  of  atoms  in  the  region  having 

velocity  in  a  small  interval  about  V,  which  is  plotted  on  the  horizontal 

axis.  Figure  10a  shows  the  distribution  functions  for  sixteen  planes 

in  the  vicinity  of  the  front  (planes  305-320),  whereas  Figure  10b  shows 

those  for  a  region  far  behind  the  front  (planes  65-80).  For  each  sample 

the  number  of  atoms  in  the  region,  n  ,  is  128.  Figure  10a  shows  clearly 

that  the  distribution  functions  near  the  front  differ  appreciably  from 
a  true  Maxwellian  curve.  Specifically,  as  a  result  of  the  planar  os¬ 
cillations  more  particles  are  found  near  extreme  values  of  the  velocity 
than  a  Maxwellian  function  predicts.  The  functions  in  Figure  10b,  on 
the  other  hand,  appear  to  be  more  nearly  Maxwellian  and  were  found,  from 
other  plots,  to  be  fairly  constant  in  time.  The  fact  that  they  are  much 
narrower  than  the  distribution  functions  of  Figure  10a  is  consistent 
with  our  observation  that  the  thermal  energy  far  behind  the  shock  front 
is  significantly  less  than  that  at  the  front.  For  comparison,  we  have 
plotted  in  Figure  11  the  velocity  distribution  function  in  the  uncom¬ 
pressed  lattice  ahead  of  the  shock.  The  Maxwellian  character  of  these 
curves  is  again  obvious  although,  of  course,  they  are  narrower  than 
those  of  Figure  10b  because  of  the  lower  temperature. 

The  above  considerations  do  not  rigorously  prove  the  existence  of 
thermal  equilibrium  behind  the  front.  However,  from  studying  many 
similar  curves,  we  are  led  to  conclude  that,  at  least  for  practical 
purposes,  the  region  of  the  lattice  far  behind  the  front  does  approach 
equilibrium.  The  resolution  of  CMD  calculations  is,  however,  limited 
enougl,  that  one  cannot  determine  precisely  where  behind  the  front 
thermal  equilibrium  is  first  re-established.  Consequently,  we  are  unable 
to  determine  conclusively  from  our  calculations  if  the  region  of  thermal 
equilibrium  lags  behind  the  shock  front,  as  claimed  by  Tsai. 
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Figure  10.  Velocity  distribution  functions.  (at  (’lanes  305-320.  (b")  Planes  65-80. 
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5.  DISCUSSION 


In  this  final  section  we  will  be  concerned  primarily  with  summa¬ 
rizing  the  results  of  our  calculations,  discussing  the  limitations  of 
the  model,  and  comparing  our  conclusions  with  those  reached  by  other 
investigators.  We  will  also  indicate  the  relevance  of  our  work  to  other 
Army-related  problems  as  well  as  suggest  possible  future  applications 
of  CMI)  techniques. 

To  re-emphasize,  we  believe  that  shock  propagation  in  the  three- 
dimensional  discrete  lattice  occurs  in  the  following  manner.  Steady 
compression  gives  rise  to  a  spectrum  of  longitudinal  solitary  waves 
which  propagate  in  the  vicinity  of  the  shock  front  and  continuously 
interact  with  the  thermal  background  of  the  lattice.  The  longitudinal 
waves  accentuate  transverse  planar  oscillations  and  can,  if  sufficiently 
strong,  produce  coupled  longitudinal-transverse  solitary  waves.  The 
existence  of  both  types  of  solitary  waves  is  made  evident  by  the  tem¬ 
perature  overshoots  which  occur  in  all  three  Cartesian  directions  for 
strong  shocks. 

The  propagation  velocity  of  the  solitary  waves  increases  with 
increasing  amplitude  and,  thus,  they  tend  to  spread  apart  as  they  propa¬ 
gate.  At  least  for  early  times  this  spreading  effect  results  in  a  grow¬ 
ing  transition  region  behind  the  front.  If  the  solitary  waves  were 
stable  in  three  dimensions,  as  they  were  in  our  one-dimensional  calcu¬ 
lations,  the  spreading  effect  would  persist  and  a  steady  profile  could 
not  exist.  The  lack  of  absolute  stability  in  three  dimensions,  however, 
leads  to  the  distinct  possibility  that  the  profile  will  become  steady 
in  a  macroscopic  sense  in  the  long-time  limit.  Unfortunately,  our  data 
are  not  sufficient  to  settle  the  question  conclusively  for  the  times 
for  which  we  ran  the  calculations,  and  the  limitations  of  the  CMD  tech¬ 
nique  preclude  longer  calculations. 

Despite  the  fact  that  considerable  information  can  be  learned  from 
CMD  shock-wave  studies,  we  should  point  out  that  caution  in  interpreting 
the  results  should  always  be  exercised.  The  most  obvious  limitation 
of  the  method  is  the  inability  to  propagate  the  shock  front  macroscopic 
distances  into  the  lattice.  This  limitation  has,  for  instance,  pre¬ 
cluded  a  more  careful  examination  of  the  question  of  steady  state  and 
the  approach  to  equilibrium  behind  the  front.  A  second  major  limitation 
of  the  method  is  the  inability  to  treat  lattices  having  reasonably 
large  cross  sections.  It  is  possible,  therefore,  that  phenomena  which 
occur  in  CMD  simulations  may  be  a  direct  consequence  of  the  small  cross 
sections  used  and  may  not  occur  in  systems  with  more  realistic  transverse 
dimensions.  For  example,  in  the  present  calculations  the  number  of 
atoms  in  a  plane  is  sufficiently  small  (32,  at  most)  that  the  average 
planar  velocity  will  not  be  zero  even  in  the  initial  thermal  background. 
Interaction  of  the  longitudinal  solitary  wave  with  these  background 
planar  oscillations  enhances  the  decay  of  the  solitary  waves  and,  in 
some  cases,  leads  to  the  generation  of  coupled  longitudinal  -  transverse 


solitary  waves.  As  the  number  of  particles  in  the  cross  section  is 
increased,  the  amplitude  of  the  planar  velocity  in  the  thermal  back¬ 
ground  is  decreased.  Therefore,  in  a  crystal  having  a  considerably 
larger  cross  section  the  longitudinal  solitary  waves  may  in  fact  be 
more  stable  than  those  observed  here,  and  the  appearance  of  coupled 
solitary  waves,  if  it  occurs  at  all,  will  be  significantly  delayed. 

In  spite  of  the  deficiencies,  we  believe  that  CMI)  calculations  provide 
a  reasonable  qualitative  picture  of  what  happens  in  real  solids  sub¬ 
jected  to  shock  compression.  In  looking  at  details,  however,  one  must 
be  exceedingly  careful  and  for  this  reason  we  are  unable  to  provide  a 
definite  answer  to  the  question  of  whether  or  not  the  shock  profile  is 
steady . 


We  would  now  like  to  turn  our  discussion  to  the  definition  of 
temperature  behind  the  shock  front  and,  hopefully,  to  shed  some  light 
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upon  a  controversy  between  Tsai  ’  and  Paskin  that  has  existed  in  the 
literature  for  some  time.  The  disagreement  arises  only  over  the  z 

component  of  the  temperature.  To  determir  this  quantity,  Paskin1* 
first  defines  the  z  component  of  the  planar  temperature  by 


(5.1) 


where  the  sum  extends  only  over  atoms  in  a  particular  plane  normal  to 

the  z  axis.  Thus,  n  is  the  number  of  atoms  in  the  plane  and  <  V  > 

p  z 

represents  the  average  instantaneous  velocity  of  the  plane  in  the  z 

direction.  The  z  component  of  the  temperature  in  a  region  is  then 

obtained  by  averaging  (Oj)  over  the  planes  in  the  region.  Presumably, 

z 

the  rationale  for  this  perhaps  unusual  definition  is  that  one  thereby 
does  not  consider  the  energy  associated  with  the  planar  oscillations 
near  the  front  as  part  of  the  random  thermal  energy  or  temperature. 
Tsai,  on  the  other  hand,  has  insisted  upon  using  the  definition 


where  the  sum  now  extends  not  only  over  atoms  in  a  given  plane,  but  also 
over  many  planes  (i.e.,  making  up  a  region  R) .  We,  for  comparison,  have 
employed  the  definition 


(5.3) 


where  <  \ ^  >  here  represents  the  actual  average  velocity  in  the  region 
under  consideration. 


The  first  point  we  wish  to  make  is  that  none  of  the  above  three 
definitions  is  "wrong".  Since  temperature  has  a  precise  meaning  only 
in  the  context  of  cqui  1  ihrium  statistical  mechanics,  one  is  free  to 
define  it  arbitrarily  in  the  nonequilibrium  region  behind  the  shock 
provided  the  definition  reduces  to  the  standard  one  in  equilibrium. 

All  three  definitions  satisfy  this  criterion  so  the  only  point  open  to 
discussion  is  which  definition  is  the  most  reasonable. 

Tsai's  definition  of  temperature  is  essentially  the  same  as  ours 
except  that  Tsai  has  assumed  that  the  translational  velocity  behind  the 
shock  is  everywhere  given  by  the  compression  velocity  1)^.  The  reason 

that  he  has  chosen  to  use  this  value  rather  than  to  calculate  the  actual 
average  velocity  in  the  region  considered  is  not  clear,  but  apparently 
he  believes  that  the  compression  velocity  is  achieved  almost  immediately 
behind  the  front.  Our  results  do  not  entirely  substant iate  this  belief. 
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We  found,  for  instance,  that  when  the  shock  was  at  the  323  plane  (see 
figure  6),  the  average  value  of  Vz  calculated  between  planes  305 

and  320  was  1.57  and  between  planes  289  and  304  was  1.83.  Both  these 
values  are  significantly  less  then  the  final  cqui 1 ibrated  value  of  2.0 

attained  still  farther  behind  the  front.  The  use  of  1J  rather  than 

P 

<  V,  >  in  the  definition  of  temperature  would  lead  to  a  higher  tempera¬ 
ture  overshoot  in  Tsai's  profiles  than  in  ours. 

1'he  principal  conceptual  difference  in  the  three  definitions,  how¬ 
ever,  is  whether  to  include  the  planar  oscillations  as  thermal  energy 
(as  Tsai  and  we  choose  to  do)  or  as  mean  flow  energy  (as  I’askin  does). 
The  choice  significantly  affects  the  temperature  profile  since  the 
overshoot  at  the  shock  front  disappears  if  the  planar  oscillations  are 
not  included  in  the  definition  of  the  temperature.  Actually,  the  issue 
is  moot  and  one  can  make  an  argument  for  either  choice.  Since  the 
planar  oscillations  arc  coherent,  well-defined  waveforms,  they  do  not, 
in  fact,  represent  the  random  motion  ordinarily  associated  with  thermal 
energy.  On  the  other  hand,  the  planar  oscillations  are  sufficiently 
narrow  that  there  is  a  substantial  difference  in  the  velocity  of  ad¬ 
jacent  planes.  The  oscillations,  therefore,  contribute  significantly  to 
the  relative  velocities  of  atoms  in  adjacent  planes  and,  in  that  sense, 
have  an  effect  similar  to  random  thermal  motion.  In  any  case,  we  feel 
that  it  is  more  important  to  understand  why  the  interplanar  oscillations 
arise  and  how  they  might  affect  the  properties  of  the  crystal  than  to 
attempt  to  categorize  them. 

A  major  objective  of  studying  shock  propagation  in  discrete  lat¬ 
tices  is  to  determine  how  the  results  compare  with  more  standard  con¬ 
tinuum-mechanical  calculations.  These  considerations  are  important  to 
a  number  of  Army-related  problems,  particularly  detonation,  as  we  have 
discussed  in  Reference  1.  Currently  used  theories  of  shock- induced 
detonation,  for  instance,  are  based  on  the  assumption  that  the  details 


of  the  shock  front  can  be  ignored.  Thus,  the  only  function  of  the  shock 
is  to  raise  the  temperature,  density  and  pressure  of  the  solid  behind 
the  shock  front  to  values  higher  than  ambient.  Chemical  reactions  are 
then  assumed  to  occur  in  the  high-temperature,  thermally  equilibrated 
region  behind  the  shock. 

In  the  present  calculations,  however,  we  have  found  that  there 
exists  a  sizable,  if  not  growing,  region  immediately  behind  the  front 
which  is  far  removed  from  thermal  equilibrium.  It  would  be  desirable 
to  determine  how  such  a  region  of  nonequilibrium  would  affect  chemical  - 
reaction  rates.  Unfortunately,  the  question  is  difficult  to  answer 
because  very  little  is  known  about  chemical  reactions  in  solids,  and 
less  still  about  how  the  rates  are  affected  by  conditions  of  non¬ 
equilibrium.  Nevertheless,  it  seems  reasonable  to  speculate  that 
reactions  are  more  likely  to  occur  in  regions  of  the  lattice  where  both 
the  thermal  and  potential  energies  of  the  atoms  are  high. 

Therefore,  to  estimate  very  roughly  the  likelihood  of  reactions, 
we  have  plotted  in  figure  12  the  sum  of  the  potential  and  thermal 
energies  behind  the  shock.  The  thermal  energy  has  been  obtained  from 
Eq.  (5.3)  except  that,  instead  of  averaging  over  16  planes  to  calculate 
0,  we  first  averaged  over  only  two  planes,  obtained  the  appropriate  6, 
and  then  averaged  the  result  over  8  consecutive  sets  of  two  planes. 

The  purpose  of  performing  the  average  in  this  manner  is  to  retain  only 
that  portion  of  the  planar  oscillations  which  leads  to  relative  motion 
of  atoms  in  adjacent  planes,  which  are  the  most  likely  to  react.  We 
did  find,  however,  that  performing  the  average  in  this  manner  led  to 
results  similar  to  those  obtained  by  simply  averaging  over  16  planes  in 
Eq.  (5.3)  except  that  the  temperature  overshoot  at  the  shock  front  was 
reduced  by  approximately  15°o.  This  relatively  small  decrease  in  the 
temperature  at  the  front  indicates  that,  at  least  for  this  case,  the 
width  of  the  solitary  waves  is  of  the  same  order  of  magnitude  as  the 
interatomic  separation.  In  fact,  it  is  easy  to  verify  from  the  planar 
trajectories  that  the  solitary  waves  near  the  front  extend  over  approxi¬ 
mately  three  lattice  planes. 

The  quantity  e'  plotted  in  figure  12  represents  the  average  thermal 
and  potential  energy  per  particle  in  regions  of  the  lattice  containing 
16  planes  along  the  z  axis.  The  temperature  was  calculated  in  the 
manner  discussed  above  and  is  related  to  the  thermal  energy  by  the 
relation 


Ejh  =  3/2  nR  0;  (5.4) 

the  zero-point  potential,  or  potential  of  the  uncompressed,  quiescent 
lattice  was  also  subtracted  from  the  definition  of  e  .  Therefore,  e ' 
was  given  by  the  relation 
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where  <t> ^ ^ ^  is  the  potential  energy  of  the  atom  (i,j,k),  viz.. 
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and  4^  is  the  total  zero-point  potential  of  the  region  under  consider¬ 
ation. 

As  can  be  seen  from  the  figure,  ey  is  fairly  constant  behind  the 
shock.  The  absence  of  an  overshoot  in  e'  results  from  the  fact  that 
the  potential  energy  just  behind  the  shock  front  is  lower  than  it  is  in 
the  equilibrated  region  far  behind  the  front  and,  in  fact,  the  difference 
is  sufficient  to  compensate  for  the  overshoot  in  the  thermal  energy. 

This  result  is  in  agreement  with  that  of  Paskin11  who  found  no  over¬ 
shoot  in  the  total  energy  density  behind  the  shock.  Therefore,  if  only 
c1  were  the  relevant  parameter  for  determining  the  probability  of 
reactions.  Figure  12  would  suggest  that  reactions  are  equally  probable 
in  the  equilibrated  and  relaxing  regions  of  the  lattice.  The  assumption 
that  equilibrium  was  attained  immediately  behind  the  front  would  perhaps 
be  satisfactory  then,  at  least  insofar  as  the  chemistry  of  the  problem 
is  concerned. 


Probably  a  more  reliable  estimate  of  the  reaction  probability, 
however,  is  the  distribution  function  of  c‘  at  various  positions  in  the 
lattice.  This  distribution  function  is  plotted  in  Figure  13  for  two 

rd 

regions  of  the  lattice  when  the  shock  front  was  at  the  323  plane 
(1^=12.80),  Figure  13a  depicts  the  function  evaluated  between  planes 

65  and  80,  presumably  a  region  considerably  removed  from  the  thermally 
relaxing  region  behind  the  shock.  Figure  13b  shows  the  function 
evaluated  between  planes  305  and  320,  immediately  behind  the  front. 

The  vertical  axis  represents  the  fractional  number  of  atoms  in  the 
region  having  energy  e/  in  the  range  shown  on  the  horizontal  axis. 

The  most  obvious  difference  in  the  two  functions  is  that  the  dis¬ 
tribution  immediately  behind  the  front  is  somewhat  broader  than  that 
obtained  in  the  equilibrated  region.  Consequently,  more  of  the  atoms 
in  the  lattice  have  extreme  values  of  ey  in  the  relaxing  region  than 
in  the  equilibrated  region.  Since  large  values  of  e1  are  more  likely 
to  lead  to  reactive  collisions  than  small  values,  one  might  speculate 
that  the  transistion  region  behind  the  front  is  more  conducive  to 
chemical  reactions  than  is  the  region  in  thermal  equilibrium. 
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The  above  analysis  is  decidedly  imprecise  and  is  intended  only  as 
crude  speculation.  We  do  not  know,  for  instance,  precisely  how  the 
absence  of  equilibrium  affects  reaction  rates  or  even  if  the  parameter 
e;  is  the  most  relevant  one  to  study.  Furthermore,  if  the  transition 
region,  which  appears  to  be  100-200  planes  at  1^=12.80,  does  not  grow 

with  time,  the  region  may  not  be  long  enough  to  affect  the  reaction 
mechanism.  At  any  rate  we  do  believe  that  our  results  cast  some  doubt 
upon  conventional  theories  of  detonation  and  suggest  a  need  for  further 
study  of  chemical  reactions  in  solids.  Of  particular  interest  is  to 
determine  how  the  reaction  rates  are  affected  by  a  nonequilibrium 
environment  and  what  the  relevant  parameters  to  study  are. 

In  future  extensions  of  the  work  undertaken  here  it  would  be 
desirable  to  make  the  model  more  sophisticated.  The  addition  of  iso¬ 
topic  impurities  would  be  a  simple  extension  which  might,  nonetheless, 
significantly  affect  the  results.  Somewhat  more  difficult  would  be  to 
treat  the  effects  of  vacancies  and  other  types  of  defects  in  the  lattice. 
Still  more  difficult  would  be  an  attempt  to  model  in  some  approximate 
way  a  molecular  crystal  in  order  to  study  the  exchange  of  energy  be¬ 
tween  intermolecular  and  intramolecular  vibrations.  The  first  of  these 
tasks  is  planned  for  the  near  future. 
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